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Abstract 

A computational strategy for the evaluation of the crystal lattice constants and cohesive energy 
of the weakly bound molecular solids is proposed. The strategy is based on the high level ab initio 
coupled-cluster determination of the pairwise additive contribution to the interaction energy. The 
zero-point-energy correction and non-additive contributions to the interaction energy are treated 
using density functional methods. The experimental crystal lattice constants of the solid benzene 
are reproduced, and the value of 480 meV/molecule is calculated for its cohesive energy. 
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I. INTRODUCTION 



In the last 50 years, weakly bound solids, including rare-gas solids and organic molecular 
solids, have been the subject of considerable research interest.-^ Many experimental mea- 
surements and theoretical predictions of the crystal structures, lattice constants, cohesive 
energies, lattice dynamics and phase behavior of weakly bound solids have been published. 
While experiments can be performed directly on the weakly bound solids, theoretical predic- 
tions are often derived from models that build upon pairwise additive interaction energies 
and include the most important refinements such as non-additive contributions and zero- 
point energy. Such information is obtained from calculations on the corresponding weakly 
bound dimers, trimers or tetramers. Therefore, a good knowledge of the intermolecular 
interactions is essential for the reliable theoretical determination of the properties of the 
weakly bound solids. 

Recently, a series of the state-of-the-art theoretical studies on rare-gas solids has been 
performed using the high-level ab initio coupled-cluster method.- 1 ^^ In these studies, the 
pairwise additive contributions to the interaction energy were represented by empirical po- 
tentials, the zero-point-energy (ZPE) corrections were calculated within the harmonic ap- 
proximation from the pairwise additive contributions, and non-additive contributions were 
determined using coupled-cluster calculations with single and double excitations and a per- 
turbative treatment of the triple excitations (CCSD(T)). 8 It was revealed that theoretical 
predictions of the structural parameters at this level of theory agree with experiment, well 
within the experimental error bars. 

It would be ideal if it were possible to apply the same treatment to weakly bound molec- 
ular solids. However, empirical potentials are not always available. Moreover, when the cor- 
responding molecule is relatively large, the high level ab initio methods, such as CCSD(T), 
are not always affordable. That is why density-functional-theory (DFT) calculations are 
often performed to determine lattice constants and cohesive energies. On the other hand, 
predicting the structural parameters of weakly bound molecular crystals using merely the 
DFT methodology is not straightforward due to the inability of the DFT to account for the 
dispersion interactions. 

In this paper, we present a combined approach to the weakly bound molecular crystals 
in which the pairwise additive contribution to the interaction energy is calculated at the 



2 



CCSD(T) level and the ZPE correction and non- additive contributions are treated purely 
at the DFT level. We illustrate our approach by calculating the cohesive energy and crystal 
lattice parameters of solid benzene, the simplest real crystalline system, in which interactions 
between aromatic molecules can be studied. In a very recent study on the solid-benzene vi- 
brational dynamics,- the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 
performed extremely well, and its results were in an excellent agreement with the inelastic 
neutron-scattering spectroscopy data. For precisely this reason, we have selected the PBE 
functional as a basic DFT tool for our study. The methodology of our approach is described 
in Section II. In Section III, the results obtained for the solid benzene are discussed, followed 
by conclusions in Section IV. 

II. METHODOLOGY 

In general, the interaction energy of a molecular crystal can be expressed as 



where the 3x3 matrix a represents lattice parameters, is a pairwise additive (two-body) 
contribution, Y^ n >3 represents the sum of many-body non-additive terms, and AE Z pe 
is a zero-point-energy correction. It has been shown that the periodic DFT methods can 
accurately describe the crystal structure (mutual orientation of molecules) and spectra of 
molecular solids. On the other hand, the DFT with local density functionals badly fails for 
binding, and therefore a higher level of theory has to be used for the calculation of cohesive 
energies and equilibrium lattice constants. Our computational methodology consists of the 
following steps: (i) a fixed-volume plane-wave DFT geometry optimization for the given 
lattice parameters a, (ii) a calculation of the two-body term using a coupled- cluster method, 
and (iii) an estimate of the many-body contributions and a calculation of the ZPE correction 
at the periodic DFT level. The pairwise additive contribution E^ represents by far the 
largest contribution to the total interaction energy, which is the reason why this term was 
evaluated at the highest level of theory feasible for solid benzene. 
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A. Plane-wave DFT calculation 

The fixed- volume structural optimizations for the solid benzene were performed with the 
lattice constants obtained from the one-parameter scaling of the experimentally determined 
unit cell (Pbca orthorhombic cell with a = 7.355 A, b = 9.371 A, and c = 6.699 A). 13 The 
calculations were carried out using the periodic plane-wave DFT with the Perdew-Burke- 
Ernzerhof (PBE) exchange-correlation functional.^ The cut-off energy of 800 eV and hard 
potentials for carbon and hydrogen atoms (ENMAX=700 eV) were employed. The Kohn- 
Sham equations were solved using a plane-wave basis set by the projector- augmented wave 
(PAW) method of Blocht^ as adopted by Kresse and Joubert.— The Brillouin-zone sampling 
was carefully checked for convergence with the number of k points. The Vienna Ab initio 
Simulation Package (VASP)^ was used for all the plane-wave PBE calculations. 

B. Calculations of the term 

The CCSD(T) calculations of the interacting benzene pairs in a crystal (for intermolecular 
distances of less than 10 A) were carried out with the augmented correlation-consistent 
valence-double-C basis set with polarization functions^ (AVDZ) using the MOLPRO 2002.6 
program suite.— The complete basis set (CBS) extrapolation was performed using a simple 
correlation-energy dependence on the basis-set cardinal number X {Ex = -Ecbs + v4X -3 ).— 
The density-fitting spin-component-scaled MP2 (SCS-MP2) method^ was employed for the 
correlation energy extrapolation (using the AVTZ and AVQZ basis sets). The Hartree-Fock 
energies were calculated using the AV5Z basis on the carbon atoms and the V5Z basis on the 
hydrogen atoms (the AV5Z* basis set). The CCSD(T)/CBS estimate E^ s f° r ^ ne benzene 
dimer was evaluated according to the formula 

pCC _ P CC 7-.MP2 , P MP2-HF , P HF 

-^CBS — -^AVDZ -^AVDZ ~>~ -^CBS ^MbZ* 

All the calculated interaction energies were corrected for the basis set superposition error 
(BSSE) using the counterpoise correction method of Boys and Bernardi.— The calcula- 
tions were performed using the benzene geometry of Gauss and Stanton^ (rcc = 1.3915 A, 
r CH = 1.0800 A). 
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FIG. 1: Asymptotic part of the intermolecular benzene-benzene potential calculated for the sand- 
wich and T-shaped structures. 

C. Asymptotic quadrupole-quadrupole and dispersion corrections to E^ 2 ' 

An asymptotic intermolecular benzene-benzene potential was derived from the coupled- 
cluster data for the sandwich (D^) and T-shaped (C2 V ) structures of the benzene dimer. 
The following functional form of the asymptotic behavior was assumed 



asympt 



(3) 



R 5 R e R s 

where R is the distance between the monomer centers of mass, Q zz is an effective (as it 
includes higher-order electrostatic interactions) quadrupole moment of the benzene molecule, 
and C*6 and Cg are isotropic parameters of the dispersion interaction. The angular factor 
uj describes an angular dependence of the quadrupole-quadrupole interaction for symmetric 
top molecules^- 1 {uj = —3 for the T-shaped structure and uj = 6 for the sandwich structure 
of the benzene dimer). One-dimensional CCSD(T)/AVDZ scans through the benzene-dimer 
potential energy surface were calculated for six points between 10 A and 20 A with a step 
of 2 A for both structures (see FigCQ). Then the parameters Q ZZ} Cq, and Cg were fitted 
to the BSSE-corrected CCSD(T)/AVDZ interaction energies of the benzene dimer with the 
following results: Q zz = -5.509(5) a.u., C 6 = -1.66(4) 10 3 a.u., and Cg = -1.7(2) 10 5 a.u. 
Eq.© was used to calculate the asymptotic correction to the solid-benzene pairwise additive 
energy for intermolecular distances larger than 10 A. 
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D. Plane-wave DFT calculation of many-body and ZPE corrections 

The sum of many-body non-additive terms ^ n>3 was approximated as the difference 
between the total crystal energy (without ZPE) and the pairwise additive contribution, both 
calculated at the the PBE/P W-800eV level. A simple cubic cell with a lattice parameter of 20 
A was employed for the benzene-pair plane-wave calculations. The calculated benzene-pair 
interaction energies were in excellent agreement with the counterpoise corrected benzene- 
dimer interaction energies obtained from the PBE/AVQZ calculations. The pairwise additive 
contributions for the intermolecular distances larger than 10 A were calculated using Eq.Q 
with Q zz = —4.85(3) a.u. and Cq = —3.87(3) 10 2 a.u. (with C$ set to zero) obtained by 
fitting the corresponding PBE/AVQZ data. 

The ZPE correction AE^pe was evaluated using zone-centered (r-point) harmonic fre- 
quencies calculated at the PBE/PW-800eV level. 

III. RESULTS AND DISCUSSION 

In Table HI the individual contributions to the cohesive energy of solid benzene are pre- 
sented, including four approximation levels of the pairwise additive contribution. By com- 
paring the first two rows, it is evident that the basis-set dependence has to be considered very 
carefully. Our CCSD(T)/CBS estimates can be directly compared with the latest theoret- 
ical benchmarks for the benzene dimer. For example, our binding energy for the T-shaped 
(C20) benzene-dimer structure overestimates the best theoretical value, calculated at the 
QCISD(T)/CBS level,— by 0.6%. The same error has to be expected for the two-body con- 
tribution in Table I, and our value is thus likely to be overestimated by 4 meV/molecule. It 
arises from Table I that when taking into account only the nearest neighbors, the two-body 
contribution is underestimated by 12% even at the CBS limit. The inclusion of the pairs up 
to 10 A explicitly (Eq. [2])) and the long-range pairs asymptotically (Eq. Ej)) was found to 
be necessary. 

The pairwise additive contribution to the cohesive energy is lowered by the many-body 
and ZPE-correction contributions thus leading to our best estimates of 480 meV/molecule 
(Ceiie) and 486 meV/molecule (CqDq) of the solid-benzene cohesive energy. The calculated 
cohesive energy compares well with the experimental values of the heat of sublimation of 
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TABLE I: Two-body, many-body, and ZPE contributions to the total cohesive energy of solid 



benzene (in meV/molecule). 

model level -E-coh 

two-body nearest neighbor CCSD(T)/AVDZ 436 

CCSD(T)/CBS 507 

two-body R < 10 A CCSD(T)/CBS 571 

two-body incl. asympt. corr. CCSD(T) 592 

many-body terms PBE/800eV -68 

ZPE corr. term a PBE/800eV -44 (-38) 

best CCSD(T)+PBE estimate 480 (486) 



°the value for CeD6 is given in parentheses 

460-560 meV/molecule.— Note that by definition the enthalpy of sublimation is equal to the 
cohesive energy only at the temperature of absolute zero. 

The sum of the many-body and ZPE- correct ion terms forms 17% of the pairwise additive 
contribution. This is only slightly higher than in the case of the heavier rare gas solids 
(Ar, Kr, Xe), where the pairwise additive contributions are lowered by many-body and ZPE 
contributions as well, and the sums of the many-body and ZPE terms are in the range of 
11-15% of the pairwise additive contributions.- 

This similarity between solid benzene and heavier rare-gas solids is also apparent when 
the role of the individual contributions to the crystal lattice constants is examined. Table [III 
contains the ratios of the theoretical lattice constants calculated at various approximation 
levels to the experimental ones for rare-gas solids and deuterated solid benzene. It is obvious 
that the role of the individual contributions is rather similar, particularly in the case of the 
heaviest rare-gas solid, solid xenon. 

The asymptotic many-body contribution to the interaction energy has been approximated 
only at the DFT level (because it was beyond our computational possibilities to perform 
a CCSD(T) calculation on the benzene trimer and tetramer). From the examination of 
the results on the deuterated solid benzene (Table HI!), ^ seems that this approximation 
was sufficient for the accurate determination of the crystal lattice constants. However, this 
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TABLE II: Two-body, many-body, and ZPE contributions to the theory /experiment ratio for the 
fee crystal lattice constant of rare-gas solids (derived from Ref.[5|) and deuterated solid benzene. 





Ne 


Ar 


Kr 


Xe 


C 6 D 6 


two-body term 


0.958 


0.982 


0.980 


0.987 


0.987 


incl. many-body terms 


0.963 


0.989 


0.987 


0.993 


0.993 


incl. ZPE corr. term 


1.001 


1.000 


0.993 


0.997 


1.000 



might also be caused by a fortuitous cancellation of the terms with different signs, because 
the asymptotic many-body long-range terms are strongly geometry- dependent, and so in 
other weakly bound molecular solids such an approximation may prove to be insufficient. 

IV. CONCLUSION 

We have presented a computational strategy for the theoretical determination of the co- 
hesive energy and crystal lattice constants of weakly bound molecular solids. This strategy 
is based on the high level ab initio (CCSD(T)) calculation of the pairwise additive contri- 
bution to the interaction energy, whereas the zero-point-energy correction and non-additive 
contributions are treated using the DFT (PBE). Having applied this strategy to solid ben- 
zene, we have been able to reproduce its experimental lattice parameters as well as predict 
its cohesive energy (480 meV/molecule for CgH 6 and 486 meV/molecule for CqD 6 ). Fur- 
ther tests and possible improvements of the proposed computational scheme are desirable, 
particularly the role of the asymptotic many-body terms in weakly bound molecular solids 
needs closer attention. More studies on the weakly bound molecular crystals are currently 
being prepared. 
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